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We describe a test particle approach based on dynamical density functional theory (DDFT) for 
studying the correlated time evolution of the particles that constitute a fluid. Our theory provides 
a means of calculating the van Hove distribution function by treating its self and distinct parts as 
the two components of a binary fluid mixture, with the 'self component having only one particle, 
the 'distinct' component consisting of all the other particles, and using DDFT to calculate the time 
evolution of the density profiles for the two components. We apply this approach to a bulk fluid 
of Brownian hard spheres and compare to results for the van Hove function and the intermediate 
scattering function from Brownian dynamics computer simulations. We find good agreement at 
low and intermediate densities using the very simple Ramakrishnan-Yussouff [Phys. Rev. B 19, 
2775 (1979)] approximation for the excess free energy functional. Since the DDFT is based on the 
equilibrium Helmholtz free energy functional, we can probe a free energy landscape that underlies 
the dynamics. Within the mean-field approximation we find that as the particle density increases, 
this landscape develops a minimum, while an exact treatment of a model confined situation shows 
that for an ergodic fluid this landscape should be monotonic. We discuss possible implications for 
slow, glassy and arrested dynamics at high densities. 

PACS numbers: 05.20.Jj, 61.20.-p, 64.70.Q- 



I. INTRODUCTION 

The structure of condensed matter is commonly 
probed by X-ray or neutron scattering techniques, that 
yield quantities such as S(k), the static structure factor, 
and its dynamical counterpart F(k,t), the intermediate 
scattering function 1 . However, in recent years, with the 
advent of modern confocal microscopes, which are able 
to characterise the structure of colloidal suspensions in 
real-space, there is an equally great emphasis on deter- 
mining the radial distribution function g(r) and its dy- 
namical counterpart G(r, £), the van Hove distribution 
function^ - — . The van Hove distribution function G(r, t) 
is a real-space dynamical correlation function for char- 
acterising the spatial and temporal distributions of pairs 
of particles in a fluid. It gives the probability of find- 
ing a particle at position r at time t, where |r| = r, 
given that one of the particles was located at the ori- 
gin at time t = 0. The intermediate scattering function 
F(k, t) is simply obtained from G(r, t) via spatial (three- 
dimensional) Fourier transform. Pair correlation func- 
tions are important because of the significant amount 
of information that they contain: Transport coefficients 
can be calculated via Kubo formulae - for example the 
diffusion coefficient D can be obtained from G(r, t) - 
and thermodynamic quantities such as the internal en- 
ergy and the pressure can be related to spatial integrals 1 
involving g(r). Whether or not a liquid is near to freez- 
ing can often be discerned from inspecting the height 
of the principal peak in S(k): it was first noticed by 
Hansen and Verlet^ that many simple liquids freeze when 
the principal peak in S(k) at k = k m , takes the value 



S(k m ) ~ 2.85. Whether a system is a glass (i.e. an amor- 
phous solid) rather than a fluid may be determined from 
the long time limit value of F(k, t), because in a fluid the 
t — » oo limit of F(k,t) is zero, whereas for a glass this 
limit takes non-zero values. This brief (and incomplete) 
survey is intended to demonstrate that both dynamical 
and static pair correlation functions are fundamental for 
characterising and understanding liquids. 

In the history of liquid state physics, fluids of hard 
spheres have proved to be an important model system 
for developing new techniques and theories. The hard 
sphere model is composed of particles interacting via the 
pair potential 

, s _ J oo r < a , . 

Vhs(r) -< [Q otherwise, 1 > 

where r is the distance between the centers of the par- 
ticles and a is the hard sphere diameter. Hard spheres 
play an important role in describing real systems, be- 
cause attractive interactions such as those present in the 
Lcnnard-Jones potential, can often be treated as a per- 
turbation to the hard sphere system^. Hence a theory 
that can successfully describe the properties of the hard 
sphere fluid forms a good candidate to work for more 
realistic systems. The hard sphere model has further 
grown in importance in recent decades due to the fact 
that Eq. fTJ) provides a good model for the effective in- 
teraction potential between colloidal particles in suspen- 
sion, in the case when the charges on the colloids are 
small or well screened - see e.g. Ref. d for an example of 
such a system. As the density of a hard sphere fluid is 
increased, the system freezes to a crystalline state, and 
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the hard sphere model has played an important role in 
developing our understanding of this phase transition. 
In contrast, although the glass transition has attracted 
much interest in recent years, it is still not completely un- 
derstood. An introduction to the vast literature on this 
subject can be found in Refs. [l|@ and references therein. 
A number of universal processes have been discovered, in- 
cluding dynamical heterogeneity^, stretched exponential 
decay of correlation functions 8 , and two-stage relaxation 
times^. In order to understand the processes involved 
in structural arrest, a number of different theoretical ap- 
proaches have been used. In particular, mode-coupling 
theory (MCT) has been successful in describing the bulk 
glass transition for hard sphere colloids-, and has been 
applied e.g. to suspension rheolog y 10 ' 11 . Nevertheless, 
alternatives to MCT have been developed 1 ^— . What is 
clear from the many studies of arrested systems, is that 
key signatures of the slow dynamics are manifest in dy- 
namical pair correlation functions. 

In our previous Rapid Communication 1 ^ a theory to 
calculate the van Hove function was proposed. For a 
bulk fluid of particles interacting via Gaussian pair po- 
tentials, comparison with Brownian dynamics (BD) com- 
puter simulation results showed that the theory is very 
reliable for determining G(r, t) for this particular model 
system. The theory is formulated for inhomogeneous 
systems and hence was also applied to investigate the 
dynamics of hard spheres confined between two parallel 
hard walls 1 ^. This approach has since been applied to 
investigate dynamics in liquid crystalline systems 1 ^. In 
the current paper we explore the theory further, and ap- 
ply it to study a bulk fluid of Brownian hard spheres. 
We present results for the self and distinct parts of the 
van Hove distribution function, G s (r,t) and Gd(r,t) re- 
spectively, and by Fourier transforming, for the inter- 
mediate scattering function F(k,t). We also display 
results for the scaled intermediate scattering function 
4>(k,t) = F s (k,t)/F s (k,t — 0) evaluated at the wave 
number ka = 2ir. This function is often the central 
object of focus of MCT—. The two stage relaxation of 
4>(k,t), that MCT predicts close to the glass transition, 
is also present in our theory. We also determine G(r, t) 
and F{k, t) using BD computer simulations and compare 
these results with those from the theory. 

Our starting point is a dynamical generalisation of Per- 
cus' test-particle approach 1 ^ for determining the radial 
distribution function g(r). Percus showed that for a fluid 
of classical particles interacting via the pair interaction 
potential v(r), that if one sets the external potential act- 
ing on the fluid it(r) = v(r), then the one body den- 
sity distribution p(r) of the fluid around the fixed 'test' 
particle is equal to the radial distribution function, mul- 
tiplied by the bulk fluid density p\ i.e. Percus showed 
that p(r) = pg(r). When using equilibrium density func- 
tional theory (DFT) 1 * 2 ^ to study a fluid, the test-particle 
method provides a useful route to obtain g(r) because 
u(r) (and hence v(r)) enters the framework explicitly. 
We also describe an alternative 'zero-dimensionality' ap- 



proach for calculating g(r). This forms a stepping stone 
in the development of the dynamical theory. 

We apply a dynamical extension of Percus' idea, 
together with dynamical density functional theory 
(DDFT)2ir^ in order to calculate the van Hove function 
G(r, t) in general inhomogeneous situations. We imple- 
ment the method using the very simple Ramakrishnan- 
Youssouff (RY) approximation for the Helmholtz free en- 
ergy functional 24 . We find that the results from the the- 
ory agree well with those from BD computer simulations 
when the fluid density pa 3 < 0.6. At higher densities 
the free energy underlying the dynamics develops a min- 
imum, corresponding to the appearance of a free energy 
barrier that must be traversed for a particle to escape 
from the cage formed by the neighbouring particles. 

In addition, we compare our results for G(r, t) to those 
obtained by assuming that G s (r, t) takes a simple Gaus- 
sian form for all times t, together with the Vineyard 
approximation 1 * 2 ^, which sets Gd(r, t) to be a simple con- 
volution of G s (r, t) and the radial distribution function 
g{r), as described in detail below. We find that in con- 
trast to the received wisdom 2 ^, the simple Vineyard ap- 
proximation is actually a fairly good approximation for 
the van Hove function for Brownian hard sphere fluids at 
low and intermediate densities. 

We also compare to an equilibrium DFT based ap- 
proach with which we are able to calculate a series of den- 
sity profiles that agree well with those from the DDFT. 
This is done by performing a constrained minimisation 
of the free energy through the judicious choice of a suit- 
able external potential to confine the test particle. This 
approach is easier to implement than the full DDFT and 
allows for the free energy landscape underlying the dy- 
namics to be mapped out and examined in detail. How- 
ever, this approach does not give any of the time infor- 
mation that the full DDFT gives - i.e. it yields the van 
Hove function with the time labels 'removed'. One of 
the advantages of this approach is that for a particular 
(parabolic) choice of external potentials, we are able to 
calculate exactly the fluid density profiles, which are pre- 
cisely those predicted by Vineyard's theory 1 ^. We dis- 
cuss the significance of this result below, after we have 
laid out the general structure of the theory and shown 
the results. 

This paper is structured as follows: In Sec. [H] we out- 
line the necessary theoretical background, including the 
definition of the van Hove function, the Vineyard approx- 
imation, DDFT, and the static test particle limit. Most 
of this section may be safely skipped by expert readers. 
In Sec. lIIII the dynamical test particle limit is introduced. 
Sec. IIVI summarises the model used and describes the 
simulation details. In Sec. [V] we describe results from the 
different dynamical approaches, the corresponding equi- 
librium approaches, and the free energy landscape. In 
Sec. IVII we make some concluding remarks. Appendix lAl 
presents an exact solution of a corresponding equilibrium 
situation. 
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II. BACKGROUND 

A. The van Hove function 

We first recall the definition of the van Hove function 
and some of its properties; for a more detailed account 
see Refs. Consider a set of N particles with time 

dependent position coordinates rj(t), where i = 1,..,N 
is the particle index, and t is time. The van Hove cor- 
relation function is defined as the probability of finding 
a particle at position r at time t, given that there was a 
particle at the origin at time t = 0: 



N N 



G(r,t) 



^(EE«(' + 'i(0)-r<(*)) 



(2) 



u=l j=i 



where (•) represents an ensemble average and S(-) is the 
three-dimensional Dirac delta function. G(r, t) can be 
naturally separated into two terms, conventionally re- 
ferred to as its "self" and "distinct" part, by discrimi- 
nating between the cases i = j and i ^ j, respectively. 
So 

= ^(E^r + r^O)-^))) 



+^(l>( r + r;(0)-r i (t)) 
G s {r,t) + G d {r,t), 



(3) 



where the self part, G s (r, t), describes the average motion 
of the particle that was initially at the origin, whereas 
the distinct part, Gd(r,t), describes the behavior of the 
remaining N — 1 particles. At t = 0, Eq. (J3]) reduces 
to the static particle-particle auto-correlation function, 
which is defined as 



1 / N 

G(r,0) = $(r) + -(£«$(r + r J (0)-r i (0)) 



N 

5(r)+pg(r), 



(4) 



where g(r) is the (static) pair distribution function. For 
the homogeneous bulk fluid p(r) = p; isotropy implies 
that the dependence is only on r — |r|. Thus, at t = 0: 



G s (r,0) = <5(r) 
G d (r,0) = pg(r). 



(5) 
(6) 



From the definitions of G s (r,t) and G d (r,t), Eq. ([3]), it 
is clear that the volume integral of these functions must 
be a conserved quantity for all times t: 



drG s (r,t) = 1, 
drG d (r,t) = N-l. 



(7) 
(8) 



The asymptotic behaviour of G(r, t) in bulk in the ther- 
modynamic limit is obtained by considering N — ¥ oo and 
volume V — > oo such that N/V = p is finite: 

lim G s {r,t) = lim G s (r,t) = 0, (9) 

r— »oo t— >oo 



lim Gd{r,t) = lim G d {r,t) = p. 



(10) 



A key quantity that we use below to characterise G s (r, t) 
is its width w(t) defined via 



(w(t)) 2 =4tt / dr r 4 G s (r,t), 



(11) 



the second moment of G s (r,t). It is often convenient 
to consider the intermediate scattering function which is 
related to the van Hove function via a spatial Fourier 
transform, 

F(k,t) = /drG(r,t)exp(-ik-r). (12) 

This quantity is directly accessible in light and neutron 
scattering experiments^. 



B. Approximating G„(r,t) 

A commonly used approximation for the self part of 
the van Hove function is to assume a Gaussian shape 1 : 



G,(r,t) 



1 



w 3 / 2 w{ty 



■ exp 



W{tf 



(13) 



where the width W(t) = J if u>(t), when w(t) is calculated 

via Eq. (fTTj) . The form (fT3|) is exact in the limits t -» 
and t — > oo for all densities when the system is fluid 1 . 
It is also exact for all times t in the low density limit 
p — > 0, where interactions between the particles can be 
neglected. There are a number of approximations for 
W(t). In molecular dynamics, at very short times t <C t c , 
where r c is the mean collision time, particles in a fluid do 
not experience collisions and therefore move freely at a 
constant velocity. This is akin to an ideal gas where the 
particle velocities follow a simple Maxwellian (Gaussian) 
distribution, giving 



W(t) = ty/2/Pm, 



(14) 



where m is the particle mass. Over longer times t t c 
the particles in the fluid undergo many collisions with 
neighbouring particles, so that the trajectory of a given 
particle is a random walk and thus its probability distri- 
bution G s (r, t) is the solution of the diffusion equation 



dG s (r,t) 
dt 



AV 2 G s (r,t), 



(15) 



where Di is the long time self diffusion coefficient. For the 
Dirac delta initial condition ([5]), the solution of Eq. (fT5|) 
corresponds to the Gaussian form (|13|) . with 



W{t) 



(16) 
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For colloidal particles, the collisions with the solvent 
atoms happen so frequently that the time scale r c is much 
smaller than all other time scales relevant for the dynam- 
ics, such as the Brownian time scale tb which is roughly 
the time for a particle to travel a distance equal to its 
own diameter, and also the typical time scale between 
collisions of pairs of colloids, r co ;. This means that we 
may set r c — ¥ and that for a low density suspension 
of colloids Eq. (IT51) holds for all times t. Thus, we may 
combine Eqs. (fT5|) and (TTB| to obtain 



G s (r,t) = (4^At)" 3/2 exp( - 



(17) 



We find below that for Brownian hard spheres this ap- 
proximation is not only reliable in the low density limit, 
but is also fairly good up to intermediate densities pa 3 < 
0.6. 



C. Vineyard approximation for Gd{r,t) 

Vineyard 2 -' suggested that one may rewrite the distinct 
part of the van Hove function as 



G d (r,t) 



J dr'g(r')H(r,r',t), 



(18) 



which is merely a redefinition of Gd{r, t) in terms of the 
unknown function H(r,r',t), which is the probability 
that if there was a particle at the origin at time t = 
and a second particle located at r', this second particle 
is later located at r at time t. Vineyard's approximation 
is to replace H(r, r', t) by G s (r — r', t), giving 



G d (r,t) = J dr'g(r')G s (r-r',t). 



(19) 



Some comments in the literature state that the Vine- 
yard approximation ignores important correlations that 
inhibit the rate at which the structure of the liquid 
breaks up and it therefore predicts too rapid decay of 
this structure 2 ^. This may indeed be the case for fluids 
undergoing molecular dynamics, but for the system with 
Brownian dynamics (over-damped stochastic equations 
of motion) that we consider here, we find that taking Eq. 
(IT71) together with Eq. (jT9"]) . is actually fairly reliable - 
in particular when the fluid is at low and intermediate 
densities pa 3 < 0.6. We will henceforth refer to Eqs. 
([TTf and (TT91 as the 'Vineyard approximation' for the 
van Hove function. 



D. DDFT and equilibrium DFT 

The dynamics of a system of N Brownian (colloidal) 
particles with positions can be modelled with the 
following set of (over-damped) stochastic equations of 
motional: 



-i <M*) 

dt 



ViU N (r" ,t) + &{t), 



(20) 



where r N = {r^; i — 1, . . . , N} is the set of particle co- 
ordinates, r^ 1 is a friction constant characterizing the 
one-body drag of the solvent on the particles, d(t) is a 
stochastic white noise term and the total inter-particle 
potential energy is 

N N 

U N (r N , t) = J2 u(r if t) + 2 E E v ^ ~ ( 21 ) 

i—l i—1 j^i 

which is composed of a one-body contribution due to 
the external potential u (which may or may not be 
time dependent), and a sum of contributions from the 
pair interactions between the particles. The time evo- 
lution of the probability density for the particle coor- 
dinates p( N )(r N ,i) is described by the Smoluchowski 
equatio n 23 i 27 : 



dt 



N 



= r^V l -[fc B TV 4 F W +V l C/ w F (Ar) ]. (22) 



The one-body density is obtained by integrating over the 
position coordinates of all but one particle: 



p{ri,t) = N J dv 2 ... J dT N P{v N ,t), 



(23) 



Integrating the Smoluchowski equation (|2"2"j) we obtain 
the key equation of DDFT2 3 .: 



d P (r,t) 
dt 



TV 



p(r,i)V 



SF[p(r,t)} 
6p(r,t) 



(24) 



where F[p] is taken to be the equilibrium total Helmholtz 
free energy functional: 



F[p(r)\ = k B T J drp(r)[ln(A 3 p(r)) - 1] 



+F cx [p(r)} + / drw(r)p(r 



(25) 



where the first term on the right hand side is the ideal- 
gas contribution to the free energy, A is the thermal de 
Broglie wavelength, F cx [p(r)] is the excess (over ideal gas) 
part of the free energy, which is in general unknown ex- 
actly, and we have suppressed the dependence on temper- 
ature T and volume V in the notation. In obtaining (|24[) 
we have made the approximation that equal-time two- 
body correlations at each time t in the non-equilibrium 
situation are the same as those of an equilibrium fluid 
with the same one-body density profile p(r,i), generated 
by an appropriate external potential 2 -^— . It has been 
shown in a variety of cases that the DDFT (|2~4"|) is reliable 
in predicting the time-evolution of p(r,t), when solved 
in conjunction with a sufficiently accurate approxima- 
tion for the equilibrium Helmholtz free energy functional 
F\ p(r)] - see for example the results presented in Refs. 

Although in the following we will not go beyond dy- 
namics that are local in time, it is worth mentioning that 
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more generally, going beyond the case of particles with 
stochastic over-damped equations of motion (|2T)j). Chan 
and FinkenS established a rigorous DDFT for classical 
fluids, showing that the time evolution of the one body 
density p(r, t) is obtained from the solution of 



dp(r,t) 
dt 

di(r,t) 
dt 



-V-j(r,i), 
= P|>(r,*)], 



(26) 
(27) 



where j(r,t) is the particle current density, and Eq. (|26[) 
represents a continuity equation for the one-body density 
p{r,t). One should, of course, expect on general grounds 
for the dynamical equations to be of this forrr^— - 
recall that Eq. (|27p is the continuity equation. However, 
the functional P[p(r, f)] that governs the time evolution 
of p(r, t) , takes a form that depends on the equations 
of motion of the particles - i.e. it depends on whether 
the particles evolve under Newtonian dynamics or have 
stochastic equations of motion such as (|20p. 

Due to the fact that in general the functional P[p(r, t)) 
in Eq. (1271) is an unknown quantity, one is prevented 
from applying our DDFT approach for calculating dy- 
namic correlation functions, and we are restricted to the 
the Brownian case (|2T)|) outlined above. The particular 
approximation used in Eqs. (|2^|) and (|27[) to obtain Eq. 

is to assume the one particle current density to be 
of the form 



j(r,t) = -rp(r,t)V 



5p(r,t) 



(28) 



Nevertheless, there is much active research aimed at go- 
ing beyond the simple overdamped case^— . 

In what follows we will relate the van Hove function to 
the time evolution of the one-body density profiles of a 
binary mixture; we therefore require the multicomponent 
generalization^ , of Eq. ((24)) : 



dpt{r,t) 
dt 



TV 



Pi(r,t)V 



SF[{Pi}] 



5pi(r,t) 



(29) 



where has the following form (cf. Eq. (|2"5]) ^a£: 

F[{ Pi \] = fc S T^ J dvp t {v)[\nk z p t {v) - 1] 

i 

+F ex [{ Pi }]+J2 f Avu,{r) Pt {Y). (30) 

i 

where the summations run over all species i. Given an 
initial set of density profiles, {pi(r, t = 0)}, we may 
employ the DDFT equations (|2T)f and (|3T)|) to calculate 
the full time evolution of the one-body density profiles 

For completeness, we recall some of the key results 
from equilibrium^^. For a given set of (one-body) ex- 
ternal potentials {iii(r)}, the unique set of one-body 
density profiles {pi(r)} are those which minimize the 



Helmholtz free energy of the system subject to 

the constraint that the average number of particles of 
each species J drpi(r) = Ni, is fixed. This is equivalent 
to an unconstrained minimization of the grand potential 
functional 



where the Lagrange multiplier pi is the chemical poten- 
tial of species i. Minimization with respect to variations 
in the density profiles yields the following set of Euler- 
Lagrange equation s 1 ' 20 : 



5F[{Pi}] = 
Spi(r) 

The Euler-Lagrange equations can be rewritten as 



Pi(r) = A 3 exp 



/3^-/3 Ul (r) + Cl (1) (r[{ Pj }]) 



where 



4 1) (r;[{ft}]) = -^^¥7¥ 1 



(32) 



(33) 



(34) 



is the one-body direct correlation functional. The set of 
density profiles that satisfy Q33[) minimize the free energy 
and are the equilibrium density profiles. When the equi- 
librium set of profiles {pi(r)} are substituted into (|3ip . 
the grand potential fl of the system is obtained. 



E. Percus' test particle limit 

Here we give a derivation of Percus' (static) test par- 
ticle limit closely following Ref. |4|| Consider a one com- 
ponent system such that the Helmholtz free energy, F, 
is given by (1231) . We are interested in the change in p(r) 
when the external potential is changed from the poten- 
tial u'(r) to the potential u(r). To this end we per- 
form a functional Taylor expansion of -F C x[p] in powers 
of Ap(r) = p(r) — p'(r). For the sake of simplicity we 
consider the change in going from u'(r) — to a spheri- 
cally symmetric external potential u(r). The variable in 
the Taylor expansion is then Ap(r) = p(r) — p b , where p b 
is the bulk density and the expansion of F cx [p] to second 
order in Ap(r) is 



F cx [p] = F cx [p b ] + / dr 



SF ex [p] 



Sp(r) 

2 J J 5p(r)6p(r>) 
-0(Ap 3 ). 



Ap(r) 



Ap(r)Ap(r') 
(35) 



Although the form of F cx is not specified, the functional 
derivatives are related to identifiable properties of the 
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system, so that evaluating them at p{r) — p b gives 



sf c 4p] 



5p(r) 
S 2 F cx [p] 



Sp(r)8p(r') 

5Q(Ap 3 ) 
Sp(r) 



= -k B Tp cx , 

= -k B Tc^(\r-r'\), (36) 
= B(r), 



where p cx is the excess chemical potential, c^ 2 '(r) is the 
(pair) direct correlation function, and B{r) is an un- 
known function that contains the higher order terms of 
the Taylor expansion. Substituting Eqs. (|33|) and (j3"6")l 
into (1251) . and then minimizing the functional with re- 
spect to variations in p(r), we obtain the following Euler- 
Lagrange equation [c.f. Eq. (|3"2"])]: 



SF 
Sp(r) 



3 J>\ 



[i = fc B Tln(AV 



(37) 



where we have separated the chemical potential p into an 
ideal-gas and an excess (over ideal-gas) contribution p ex . 
In the case of a spherically symmetric external potential, 
Eq. (|37[) may be rewritten as: 



p{r) 



exp 



(3u(r) + / dr'c(|r - r'\)Ap(r') 



-B(r) 



(38) 



For the same one-component system the bulk Ornstein- 
Zernike (OZ) equation for the total correlation function, 
h(r) = g(r) — 1, reads as follows: 



h(r) = c(r) + p b / dr'h(r')c(\r - r'\). 



(39) 



It can be shown through diagrammatic-methods^ that 
the OZ equation has the general solution 



h(r) = c(r) + ln(p(r)) + /3v(r) + b(r) 



(40) 



where v(r) is the inter-particle pair potential, and b(r) 
is the bridge function composed of the sum of all the so- 
called 'bridge' diagrams 1 -. Substituting (j4"0"|) into we 
obtain: 

g(r) = exp \-/3v(r)+ [ dr' c(\r-r'\)p b h(r')+b(r)] . (41) 



If we compare Eqs. (|38[) and (|4Tj) we find they have the 
same structure, and that they may be formally identified. 
If we set u(r) = v(r) in Eq. (|38p . it can be shown through 
diagrammatic methods 1 ' 49 that b(r) — B(r) and that 



g(r) = p(r)/p», 



(42) 



or alternatively, ph{r) — Ap(r). Thus when u(r) = v(r), 
Eqs. (l38|) and (|4Tj) become identical. Therefore we note 
that not only can the OZ relationship be derived from 



the free energy functional^, but that the equilibrium 
one-body density profile in the presence of an external 
potential u(r) = v(r) is related to the (two-body) ra- 
dial distribution function via Eq. (l4"2"j) . We should recall 
at this point that although many formal statements can 
be made about the bridge function b(r), in practice it is 
an unknown function, and all theories for g(r) constitute 
some form of approximation for 6(r)i. For example, if 
we set b(r) = B(r) =0 then is equivalent to using 
the hyper-netted chain (HNC) approximation 1 to the OZ 
equation. Furthermore, Percus^ showed that by Taylor 
expanding with different functions of p(r) one may re- 
trieve the Percus-Yevick and other closures to the OZ 
equations. This result may also be generalized to inho- 
mogeneous systems 4 ^. 



F. Zero-dimensionality route to g(r) 

We present an alternative method for calculating g(r), 
although its basis is the same key idea that underpins 
Percus' test particle limit described above: that g(r) can 
be obtained from the density profile of a fluid around a 
fixed test particle. The key difference is that instead of 
treating the test particle as a fixed external potential, 
in the zero-dimensionality route we treat the test parti- 
cle via its density distribution. The density profile of a 
particle fixed at a point (i.e. in zero-dimensional space - 
hence our choice of name for this limit) takes the form of 
a Dirac delta function. Having fixed this contribution to 
the density distribution [c.f. Eq. Q], one can then cal- 
culate the density distribution of the remaining particles 
in the presence of the test particle. Specifically, we can 
write the grand potential functional as: 

n*[pg(r)}=F id [pg(r)} + F cx [S(v) + pg(r))]-p f dvpg(r), 



(43) 

where pg(r) is the density distribution of the remaining 
particles - the quantity we wish to calculate. Note that 
here, and in what follows, p is the bulk density. The ideal 
gas term Fid does not contain the Dirac delta contribu- 
tion - we have crossed over to a system with N — 1 par- 
ticles. Since the bulk fluid density p is necessarily fixed, 
we must simply minimise 51* with respect to variations 
in g(r), giving the following Euler-Lagrange equation to 
be solved for g(r): 



5SF 
Sg(r) 



0. 



(44) 



An alternative means of calculating g(r) is to treat the 
system as a binary mixture. The test particle (which we 
label 's'), with density distribution p s (r) — S(r), is one 
species and then we use the DFT for a binary mixture 
to calculate the density profile of the remain particles 
Pd(r) in the presence of the density profile p s (r) for the 
fixed particle, treating the remaining particles as a second 
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species 'd' in the mixture. Pd(r) is the solution of 



Spd{r) 



= 0, 



(45) 



where Vfi is a modified version of Eq. (pH) where p s (r) — 
<5(r) is fixed and therefore the ideal Helmholtz free energy, 
Fid[pd], does not depend on p s (r), c.f. Eq. (|4"5|) . 

When using an approximate free energy functional, 
there is a difference between the the zero-dimensional 
limit and Percus' limit for calculating g(r). This is be- 
cause in the zero-dimensionality limit, in contrast to Per- 
cus' method, the test particle at the origin does not inter- 
act with the fluid via an external potential u(r) = v(r), 
that is identical to the pair potential, but rather via an 
approximation u*(r) to the pair potential, generated by 
the approximate density functional. We calculate be- 
low in Sec. IIV CI an explicit expression for u*(r) in the 
case of hard-spheres treated using the RY approxima- 
tion for the free energy. We will also discuss further the 
relation between Percus' test particle limit and the zero- 
dimensionality limit. 



III. DYNAMIC TEST PARTICLE LIMIT 
A. Definition 

We next extend the static test particle limit and con- 
sider the dynamical situation which allows us to calculate 
the van Hove function G(r,t). The key is the following 
observation: Consider a fixed test particle of species 's' 
(self) located at the origin; due to Percus we know that 
in this situation the density distribution of the remaining 
particles p(r) = pg(r). Now consider releasing the test 
particle at time t = and allowing it to move through 
the fluid. When this happens its probability (density) 
distribution p s (r, t) changes from a Dirac delta function 
(at t = 0) to a distribution with a non-zero value for some 
points away from the origin. If we now recall the defini- 
tion of the function G s (r, t) in Eq. ([3]) , we see that G s (r, t) 
gives the probability that a particle initially located at 
the origin has moved a distance r away from the origin af- 
ter time t. Therefore, in this situation, p s (r,t) = G s (r,t) 
for all times t > 0. Similarly, if we consider how the 
remaining particles redistribute themselves as the test 
particle moves away from the origin, we see from Eq. ([3]) 
that the probability of finding any one of these parti- 
cles a distance r from the origin at time t is given by 
Gd{r, t). We label the remaining particles as being parti- 
cles of species 'cf (distinct) and having the density profile 
Pd(r,t). Thus, as in the static test particle case, we may 
connect the two parts of the van Hove function with the 
density profiles of a two-component system of species s 
and d: 



G s (r,t) 
G d {r,t) 



Ps(r,t), 
Pd(r,t), 



where species s is composed of only one particle, the test 
particle, and J dr p s (r, t) = 1, so that Eq. ((7]) is satisfied. 
We may therefore formally set the pair potential for inter- 
actions between species s particles v ss (r) = 0. The den- 
sity profile for species d must satisfy the normalization 
constraint ©, and the self-distinct and distinct-distinct 
pair potentials must be identical, v s d{r) — Vd d {r) — v(r). 
This is equivalent to modelling a one-component system, 
but treating one particle separately from the rest. Recall 
that at time t = we know the test particle's position 
exactly from Eq. ([5]) and combining this with Eqs. ([6]), 
(H2) and flU we obtain 



G s (r,t 
G d {r,t 



0) 
0) 



Ps(r,t = 0) 
Pd{r, t = 0) 



S(r), 
pg(r)- 



(47) 



The connections made in Eq. (|4l)|) between the self and 
distinct parts of the van Hove function and the density 
profiles p s (r, t) and Pd(r, t) in the dynamical test particle 
limit described above are conceptually important. How- 
ever, we have merely shifted the problem of how to de- 
termine G(r, t) onto the problem of how to determine the 
time evolution of the two coupled density profiles p s (r, t) 
and p d (r, t). The solution that we use in this paper is to 
use DDFT, i.e. equations (f!?!)|) are integrated forward in 
time with Eqs. (|47[) providing the initial time, t — 0, den- 
sity profiles. The resulting time series of density profiles 
gives the self and distinct parts of the van Hove function 
through Eq. (|46D . Henceforth, we refer to this as the 
'dynamical test particle' theory. 



B. Approximate solution 

Before proceeding to the results that we obtain from 
following the calculation scheme described above, it is 
worth examining an approximate analytical solution that 
may be obtained as follows: From Eqs. ([Ml), (© and (jM)) 
we may write the DDFT equations for the two density 
profiles p s (r,t) and p d (r,t) as: 



d Pi (r,t) 
dt 



DV 2 Pl (r,t) + TV- L(r,t)V Cl (1) (r,t)l , (48) 



where i = s,d and the diffusion coefficient D — ksTT. 
If we set the second term on the right hand side to zero 
and we set D = Di then we obtain Eq. (fTsj) for p s (r, t) = 
G s (r, t) and thus the solution to the DDFT for the self 
part of the van Hove function in this limit is the Gaussian 
form in Eq. (|17p . Similarly, for species d, when we assume 
that the second term on the right hand side of Eq. 
can be neglected, then we obtain: 



dpd(r,t) 
dt 



DV 2 Pd (v,t). 



If we now assume the Vineyard approximation 



(46) 



Pd(r,t) = J dr'g(r')p s 



(\r-r'\,t) 



(49) 



(50) 
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[c.f. Eq. (|19p], then after Fourier transforming Eq. (14T)1) . 
together with Eq. (1501) we obtain 

= -k 2 Dg(k)Pd(k,t), (51) 

where g(fc) is the Fourier transform of the radial distri- 
bution function g(r) and /3d(fc, t) is the Fourier transform 
of p d (r,t). Dividing both sides of Eq. (|5i~j) by and 
then taking the inverse Fourier transform we obtain Eq. 
(fTS")) for p s (r,t) = G s (r,t). Thus the Gaussian form in 
Eq. (fl7|) together with the Vineyard approximation ((50)) 
for pd(r, t), together form a self consistent solution to the 
DDFT equations in the dynamical test particle limit, in 
the case where we can neglect the contribution from the 
second term on the right hand side of Eq. (|4"5)) . This 
term is zero in the ideal-gas limit when the excess contri- 
bution to the free energy F cx — 0, in Eq. (|3U| . or when 
p — > 0. However, we find below for hard spheres that 
this approximation is reliable well beyond the ideal gas 
regime, which suggests that in the test particle limit ci 1 ^ 
and in Eq. (|48|) must both be slowly varying (almost 
constant) functions, so that their gradients are small. 

IV. MODEL FOR HARD SPHERES 

A. Simulation method 

In order to provide benchmark results, we calculate the 
van Hove function by integrating the equations of motion 
([2"U| using standard BD computer simulations™. In order 
to apply the algorithm we model the hard spheres with 
a steep continuous pair potential: 

pv(r) — < X . (52) 

w 10 otherwise. v ' 

We solve Eqs. (|2"0"]) using the Euler forward algorithm 
using a time step 5t = 1 x 10~ 5 tb; recall that the 
Brownian time tb = <J 2 /D, where D = TksT is the 
Stokes-Einstein diffusion coefficient. The random forces 
Q in Eq. (12U1) mimic the interaction between particles 
and solvent, and are sampled from a Gaussian distribu- 
tion with zero mean and variance 2DSt. The simula- 
tions are carried out using iV=1728 particles at densities 
pa 3 = N(a/L) 3 = 0.2,0.4,0.6,0.8, and 1 in a cubic box 
of volume L 3 . 

After an equilibration time of 50 tb, we sampled the 
distribution functions G s (r,t) and G d (r,t) at the times 
£/tb=0.01, 0.1, and 1. The distribution functions are 
averaged over all possible time intervals t/rs of a single 
simulation run. The total simulated times are 2tb, 20t£, 
and 200tb for the short, medium and long time intervals, 
respectively. The scaled intermediate scattering function 
is calculated from the density autocorrelation function 
in Fourier space, 4>{k,t) = (n fc (t)n_ fc (0))/(n fc (0)n_ fe (0)), 
where nk{t) — YliLi exp(— ik • rj(t)) are the Fourier com- 
ponents of the local number density. 



B. The excess free energy functional 

In order to implement the dynamical test particle limit 
we must (as is almost always the case in density func- 
tional theory calculations) select an approximation for 
the excess part of the Helmholtz free energy functional, 
Eq. (f5U|) . We use the RY functional^, which is obtained 
from the two-component generalisation of Eqs. (1351) and 
(I5B1 by neglecting all terms of order 0(Ap 3 ) and higher. 
We obtain: 

F cx [p s ,Pd\ = Vf eK (p) 

+ f'Jj>) { / dr( Pd (r) -p)+ J dr Ps (r)| 

-IJdrj dr'c(\r-v'\)Up d (r)-p)(p d (r')-p) 

+ (p d (r) - p) Ps (r') + p s (r)(p d (r') - p) J (53) 

where / ox (p) is the bulk excess free energy per unit vol- 
ume, V is the volume of the system, f^ x = df cx /dp and 
c(r) is the bulk pair direct correlation function of the 
hard sphere fluid with bulk density p. We use / ex and 
c(r) as given by the Percus-Yevick approximation 1 . 

Eq. (|53l) is perhaps the simplest functional that one 
could use to calculate hard sphere fluid density profiles. 
Our reasons for using this functional are: (i) The struc- 
ture of the functional is relatively simple (and as a con- 
sequence is widely used within liquid state approaches). 

(ii) Within the RY functional it is straightforward to ne- 
glect the species s intra-species interactions which is nec- 
essary to ensure that p s (r,t) represents a single particle. 

(iii) Finally, the RY functional was the first functional to 
correctly reproduce freezing phenomena in hard spheres. 

C. Static structure of the fluid 

In Fig. [1] we display the radial distribution function and 
static structure factor for a bulk fluid of hard spheres for 
the densities pa 3 = 0.2, 0.4, 0.6, 0.8 and 1. We show 
results obtained from BD simulations, together with the 
results from the static test particle limit using the RY 
approximation for the excess free energy (1531) . One can 
observe that for low and intermediate densities pa 3 < 0.6, 
the test particle results are in good agreement with those 
from the simulations. However, as the density is in- 
creased, the test particle results become less reliable. We 
see in Fig. [IJb) that the theory overestimates the contact 
value, g(r = cr + ). This in turn leads to the discrepan- 
cies in the static structure factor at high densities in Fig. 
[ljd); S(k) is obtained by Fourier transforming the data 
in (b). The overall conclusion to be drawn from Fig. [1] is 
that the test particle method combined with the RY func- 
tional provides a reliable description of the fluid structure 
for low and intermediate densities, but at higher densities 
pa 3 > 0.6, the theory is only qualitatively correct. 
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FIG. 1: (a) and (b) display radial distribution functions, and 
(c) and (d) static structure factors for a bulk fluid of hard 
spheres at densities pa 3 — 0.2, 0.4, 0.6, 0.8 and 1. Parts (a) 
and (c) are obtained from BD simulations and (b) and (d) via 
Percus' static test particle limit using the RY functional. 



We return now to the discussion of the somewhat sub- 
tle issues concerning the relation between Percus' test 
particle limit and the zero-dimensionality limit. The re- 
sults from these two calculations are not the same when 
one uses an approximate expression for the free energy, 
such as the RY functional (|55)l . Combining Eqs. ([55)1 and 
(|33|) . we obtain the following: 



p d {r) = A 3 exp (3fi d - PfiJjp) 
+ J dr'c(\r-r'\) Ps (r') 
+ [ dr'c(|r-r'|)(p d (rO-p) 



(54) 



Now, recall that a test-particle calculation involves fixing 
one of the particles at the origin, treating it as an exter- 
nal potential, and then determining the density profile of 
the fluid (species d) under the influence of this external 
potential. Doing this, using the RY approximation for 



F ex [p s ,Pd}, Eq. ([53]), we obtain: 

p d (r) = A -3 exp (3pi d - Pf^(p) - (3u{r) 
+ ( dv'c{\v-v'\){p d {r')-p) 



(55) 



Comparing equations ((54)) and (j55|) , we see that Eq. ((54]) 
is merely Eq. (]5"5l) with the external potential j3u(r) — 
j3v(r) replaced by the effective potential: 



dr' C (|r-r'|)p s (r'). 



(56) 



In the w = limit, when p s (r) — S(r), we then obtain: 

Pu*(r) = -c(r). (57) 

One consequence of this random phase-like approxima- 
tion is that the core condition is violated. The degree to 
which the core condition is violated could be used as an 
indicator towards the reliability of any approximate free 
energy functional. 

In the remainder of this paper we will display results 
and distribution functions that are derived from Percus' 
test particle results used as initial condition, though we 
will draw attention to results from the zero dimensional- 
ity route where appropriate. 

V. RESULTS 
A. Dynamic approaches 

In Fig. [2] we display the two parts of the van Hove 
functions, G s (r, t) and G d (r, t), measured in the BD sim- 
ulations for fluid densities pa 3 = 0.2, 0.4, 0.6, 0.8 and 1. 
The different curves correspond to the times t/rs = 0.01, 
0.1, and 1. G s (r, t) appears to have a near-Gaussian form 
for all times and densities. For short times G d (r,t) ex- 
hibits a correlation hole for r < a and it is highly struc- 
tured for (larger) r > a. At later times the structure 
in G d (r, t) diminishes and the correlation hole becomes 
'filled in'. Recall that Eq. ((Ill]) defines the long time limit. 
Increasing the density beyond pa 3 = 1, we find that the 
simulated system crystallises onto a regular lattice and 
there is no evidence of glass forming behaviour. 

In Fig. [3J we display the one-body density profiles, 
p s (r,t) and p d (r,t), from the DDFT dynamical test par- 
ticle method for bulk fluid densities pa 3 = 0.2, 0.4, 0.6 
and 0.8. As initial condition, p d (r,t = 0) = g(r), we 
have used Percus' test particle method for calculating 
g(r), as shown in Fig. Q] The results in Fig. |3J corre- 
spond to the same times as the BD curves displayed in 
Fig. [2] namely t/rs = 0.01, 0.1, and 1. Comparing the 
BD results in Fig. d with the DDFT results in Fig.[3Jwe 
observe that for densities per 3 = 0.2, 0.4, and 0.6, there is 
good qualitative agreement between the simulation and 
DDFT results. The p d {r, t) results show a similar amount 



10 




Ha rla 



FIG. 2: The self and distinct parts of the van Hove distribu- 
tion function, G 3 (r,t) and Gd{r,t), for a hard sphere fluid, 
measured in BD simulations at densities: (a) per 3 — 0.2 (b) 
per 3 = 0.4 (c) per 3 = 0.6 (d) pa 3 = 0.8 (e) pa 3 = 1. The results 
are plotted for times i/rs=0.01 (solid line), 0.1 (dashed line) 
and 1 (dotted line). In the semi-logarithmic scale of G s (r,t) 
versus r a Gaussian appears as a parabola. The Gd{r,t) re- 
sults are shown on a linear scale. 



of structure as the Gd(r, t) results, and p s (r, t) has a very 
similar magnitude and range as G s (r, t) , although for 
per 3 = 0.6, p s (r,t) shows some departures from the al- 
most Gaussian shape observed in the simulation results, 
particularly at t/rs = 1. For per 3 = 0.8 we find that the 
dynamic test particle method predicts that after a short 
time t/rs ~ 0.1 the density profiles p s (r,t) and pd{r,t) 
cease to change with time and that the system becomes 
'arrested'. One could interpret this state as the tagged 
particle remaining localised within the cage formed by 
the neighbouring fluid particles. We discuss the signifi- 




FIG. 3: The V and 'd' density profiles, p s (r,t) and pd(r,t), 
obtained from the dynamical test particle theory, for densi- 
ties: (a) per 3 = 0.2 (b) pa 3 = 0.4 (c) per 3 = 0.6 (d) per 3 = 0.8. 
The results are plotted for times t/rs =0.01 (solid line), 0.1 
(dashed line) and 1 (dotted line). In (d), after a short time, 
the system reaches an 'arrested state', where the density pro- 
files no longer evolve in time and the width of p s (r,t — > oo) 
remains finite. 

cance of this phenomenon in Sec. IVII 

In Fig. [4] we show the van Hove functions calculated 
using the Vineyard approximation, Eqs. (TT71) and (fTI)|) 
with Di = D = k B TT, for fluid densities pa 3 = 0.2, 
0.4, 0.6, 0.8 and 1 and times t/r B = 0.01, 0.1, and 1. 
As in the DDFT we have used g(r) calculated using the 
RY functional and Percus' test particle method, although 
one could use g{r) obtained from any reasonable method, 
including g{r) measured in the BD simulations. Compar- 
ing the Vineyard results to the BD simulation results in 
Fig. [21 we find that there is reasonably good agreement 
between the two. The form of G s (r,t) is fixed to be 
Gaussian, so there is good agreement with G s (r, t) from 
the BD simulations, though it is clear that the width of 
G s (r, t) increases more rapidly in the Vineyard approxi- 
mation. For densities per 3 < 0.8 there is a similar amount 
of structure present in Gd(r, t) for r > a in the Vineyard 
approximation as in the simulation results. However, for 
per 3 = 1 the Vineyard approximation does not exhibit the 
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FIG. 4: The self and distinct parts of the van Hove distribu- 
tion function, G 3 (r, t) and Gd(r, t), calculated using the Vine- 
yard approximation for densities: (a) pa 3 = 0.2 (b) pa 3 = 0.4 
(c) pa 3 = 0.6 (d) pa 3 = 0.8 (e) pa 3 = 1. The results are 
plotted for times t/rs=0.01 (solid line), 0.1 (dashed line) and 
1 (dotted line). 



same degree of structure that is present in the simulation 
results. 

In Fig. [S]we compare the width, w(t), of the self part 
of the van Hove function, G s (r, t) , obtained from (i) BD 
simulation results, (ii) dynamical test particle limit, and 
(iii) the Vineyard approximation. In the Vineyard ap- 
proximation the time dependence of w(t) is defined by 
Eq. (|16l) and does not depend on density, so there is only 
one master curve. This is because in the same way as in 
the dynamical test particle limit, we set Di = D, where 
D = ksTT is the short time diffusion coefficient, which 
is strictly only equal to the long time self diffusion coef- 
ficient Di in the limit p — > 0. Since w(t) cx ^/t, on the 
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FIG. 5: The width w of the self part of the van Hove function, 
defined in Eq. ((TTJ , measured in the BD simulations for pa 3 = 
0.2 (+), 0.4 (x), 0.6 (*), 0.8 (□), and 1 (o), and from the 
dynamical test particle theory, for pa 3 = 0.2 (long dashed 
line), 0.4 (short dashed line), 0.6 (dashed-dotted line), 0.7 
(solid line), and 0.8 (dotted line). Also shown is the result 
from the Vineyard approximation (thick solid grey line). In 
the main panel both axes are logarithmic; the inset displays 
the same results on linear scales. 



double logarithmic scale in Fig. [5] this is represented by 
a straight line with gradient 1/2. We find that the sim- 
ulation results are also approximately linear in this rep- 
resentation for all densities considered, but that there is 
a slowing down effect as density is increased, due to the 
fact that Di decreases as the fluid density is increased 
and is no longer equal to D. For pa 3 — 0.2, the simula- 
tion results are close to the Vineyard result. As the bulk 
density is increased, the BD results move away from this 
line. 

The dynamical test particle results for w(t) in Fig. \5\ 
exhibit a much stronger dependence on density. At low 
densities the curves are similar to the Vineyard and the 
BD results, but as the density is increased the w(t) curves 
show a slowing down, and then (unphysical) speeding 
up of the dynamics, unlike that seen in the BD results. 
This slowing down is greatly exaggerated so that the 
DDFT curve for pa 3 = 0.6 is similar to the BD result at 
pa 3 = 0.8. Furthermore, the DDFT curves for pa 3 = 0.7, 
and 0.8 have no counterpart in the simulation results. 
For pa 3 — 0.7 the w(t) curve shows extremely exagger- 
ated slowing down and speeding up. We believe that the 
unphysical speeding up for t/Ts > 10 1 is due to the fact 
that the DDFT incorrectly sets the long-time diffusion 
coefficient Di equal to the short time diffusion coefficient 
D, so that as the particle escapes the cage of neighboring 
particles, it is forced to "catch-up" to give the incorrect 
long time behavior. Note that from the Smoluchowski 
equation (22) it can be shown^i that w(t) must be sub- 
diffusive for intermediate times, and that the long time 
diffusion coefficient must be smaller than the short time 
one, a feature which is well established in Brownian dy- 
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FIG. 6: Intermediate scattering function F(k, t) as a function 
of the scaled wave vector ka, obtained by a spatial Fourier 
transform, Eq. (|12|) . of the BD simulation results for the van 
Hove functions displayed in Fig.fJ 



namics simulations and experiments^— . The curve for 
pa 3 = 0.8 shows that the system slows down so much 
that the dynamics are arrested, so that w(t oo) is fi- 
nite, as one would infer from the density profiles shown in 
Fig. [31 We postpone a discussion of the possible physical 
implications to Sec. IVII 



For completeness, we also plot the intermediate scat- 
tering function F(k,t). In Fig. [5] we display results from 
BD computer simulations, and in Fig. [7] the results from 
the DDFT. We find that for pa 3 = 0.2, 0.4 and 0.6, the 
results from both approaches exhibit very similar struc- 
ture. However, since in the DDFT the dynamics become 
arrested at pa 3 — 0.8, so F(k,t) becomes arrested after 
a very short time, unlike the BD simulations result. In 




FIG. 7: Intermediate scattering functions F(k, t), obtained by 
a spatial Fourier transform, Eq. (|12[1 . of the dynamic test par- 
ticle density profiles, p s (r,t) and pd(r,t), displayed in Fig. [3] 



Fig. [8] we plot the scaled intermediate scattering function, 



-KM) 



F s (k,t) 
F s {k,t = 0)' 



(58) 



for fixed ka — 2ir, obtained from the Vineyard approxi- 
mation and compare to the BD simulation results (Fig. 
Ela)) and the DDFT results (Fig.Eft))). At the lower den- 
sities the BD simulation results and the DDFT results are 
close to the Vineyard approximation and both show some 
slowing down with density. At the higher densities the 
BD results continue to show a steady decay. However, in 
the DDFT results for pa 3 = 0.7 we see 4>(ka = 2tt, t) de- 
cays in two stages over a much longer time. For pa 3 = 0.8 
the arrested dynamics cause <fi(k/j = 2it, t) to remain fi- 
nite in the limit t — > oo. 



B. Relating dynamic to static density profiles 

One may connect DDFT and equilibrium DFT by find- 
ing the unique set of effective external potentials {ui(r)} 
that in equilibrium generate the same set of density pro- 
files as obtained in the dynamic approach at a particular 
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FIG. 8: Scaled intermediate scattering function <f>(ka = 2ir, t) 
as a function of time t/Vs, calculated from the Vineyard ap- 
proximation (thick gray line), compared to (a) BD simulation 
results, and (b) the dynamical test particle method. The 
dynamical test particle results exhibit slowing and arrested 
dynamics for per 3 = 0.7 and 0.8. 

time t. These potentials represent the net effect of fi- 
nite time and limited diffusion preventing the fluid from 
finding the structure that minimises the system free en- 
ergy. For the two-component system considered here, the 
two external potentials fiu^r, t) may be recovered, up to 
an overall time-dependent additive constant ftpi(r,t), by 
rearranging Eq. ([53]) : 

Pm(r,t) - pm =4 1) (r;[p s ,p d \)-ln[A 3 p i (r,t)}, (59) 

where p s (r,t) and pd{r,t) are the solution of the DDFT 
at time t. In Fig. [9] we plot these external potentials cor- 
responding to the density profiles from DDFT displayed 
in Fig. [3J We find that at the lowest densities per 3 = 0.2 
and 0.4, the shape of u s {r,t) is approximately parabolic 
for all times and distances r. As the fluid density is in- 
creased, u s (r,t) departs from the parabolic form. For 
pa 3 =0.8 the curves are still parabolic at large r, but 
at small r they become distorted. We find that Ud{r 7 t) 
does not vary significantly with density. At short times 
it is dominated by strong repulsion within the hard-core 



FIG. 9: External potentials, u s (r, t) and Ud(r,t), required 
to generate equilibrium density profiles, p s (r,i) and pd(r,t), 
identical to those obtained from the dynamical test particle 
approach - see Fig. [3] The external potentials also include 
an overall time-dependent additive constant which is not in- 
dicated. Note the logarithmic x-axis for u s (r, t). At low den- 
sities u s (r, t) is approximately parabolic, as indicated by the 
straight line (dashed-dotted). At high densities u s (r, t) is dis- 
torted at small r, but still parabolic at large r. 



diameter, r < a. Recall that in order to calculate g(r), 
which corresponds to t = 0, we chose to use Percus' test 
particle method, and hence have introduced an external 
potential equal to the hard sphere potential. The strong 
repulsion found for short times is a remnant of this exter- 
nal potential. As t increases, the strength of this repul- 
sion decreases and becomes almost zero for t/rs = 1— ■ 



C. Corresponding equilibrium approach 

Having established the form of the external poten- 
tials necessary to create equilibrium fluid density profiles 
equal to the profiles calculated using the dynamical test 
particle method, we now seek a simple approximation 
for this set of external potentials, to allow us to easily 
calculate equilibrium density profiles that mimic the dy- 
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namic profiles. In other words, we seek to determine the 
full van Hove function when G s (r, t) has a given width 
w, without calculating the entire preceding time series of 
profiles. In doing this we lose time t as a function ar- 
gument and instead we must 'label' the density profiles 
with w. In what follows, we will disregard the associated 
problem of relating w to time t. 

We parametrise the external potentials using a sim- 
ple functional form. Firstly, we assume that u s (r,w) is 
parabolic for all widths: 



pu s (r,\) = Ar 2 



(60) 



where A is the strength of the confining potential, and 
w is now an unknown function of A. For the external 
potential that acts on species d we consider two options. 
The first is to assume that Ud(r, w) = 0. In this case, it is 
possible to solve exactly for the equilibrium distribution 
functions and the free energy, as outlined in Appendix El 
We find that the species s density profile is a Gaussian, 



Pa(r, A) 



exp(— f3\r 2 

(7r/^A) 3 /2 



(61) 



where the dependence of the width w on A is, 

w(X) = (2A/3)" 1 / 2 , (62) 

and the d profile is given by a convolution of the radial 
distribution function g(r), together with the Gaussian 
profile p s (r), and multiplied by p: 



Pd (r)=P dr'p s (r')ff(|r-r'|). 



(63) 



These distribution functions are identical to those from 
the (dynamic) Vineyard approach. 

The second approach that we consider is to calculate 
the density profiles without defining the external poten- 
tial Ud(r 7 w) at the outset of the calculation. Instead, we 
determine this potential self consistently 'on-the-fly' as 
part of our iterative numerical solution routine, based on 
the following considerations: Firstly, recall the normali- 
sation constraints on the van Hove function in Eqs. ([7]) 
and @. In order to satisfy the normalisation constraint 
(JZJ) on the density profile for the single tagged s particle, 
we introduce a Lagrange multiplier p, s . One may also 
consider A to be a Lagrange multiplier that enforces the 
width constraint on the profile p s (r). In our calcula- 
tions the value of p s is determined on-the-fly by enforcing 
Eq. at each step of our iterative routine. However, 
we are not able to do the same for the density profile of 
the remaining d particles, because we also must have 



j(r, w) — >■ p, as r 



oo. 



(64) 



Multiplying pd(r) by a single factor breaks this condition, 
so we are not able to simply enforce (|8) at each step 
of our iterative routine in the same way as we do for 
p s (r). The condition in Eq. (|64p implies that we require 




FIG. 10: The density profiles p s (r,w) and pd{r,w), calculated 
using the equilibrium DFT. The curves are calculated at the 
densities: (a) pa 3 = 0.2, (b) per 3 = 0.4, (c) per 3 = 0.6 and 
(d) per 3 = 0.8. The curves are chosen so that the widths w 
of p s (r,w) correspond to the same widths of G s (r,t) at times 
£/tb=0.01 (solid line), 0.1 (dashed line) and 1 (dotted line) 
displayed in Fig. [2] 



an a priori unknown inhomogeneous external potential, 
Ud(r,w), with the property that Ud(r,w) — > as r —> oo. 
This may be achieved by scaling the quantity pd(r) — 
p (instead of scaling pd(r) itself) at each step, so that 
Pdir) satisfies both (JS| and (|64p. Once convergence of the 
numerical procedure is achieved, one may then inspect 
the effective external potential Ud(r,w) by substituting 
the resulting density profiles into Eq. (j5T)j) . 

The density profiles calculated using this equilibrium 
method are shown in Fig. 1101 where we plot pd(r,w) and 
p s (r, w) having widths identical to those of the van Hove 
functions from simulation, displayed in Fig. [5] Note that 
we consider only the densities, pa 3 = 0.2, 0.4, 0.6 and 0.8. 
These equilibrium profiles have been calculated using a 
normalisation constant taken from the approximation for 
g{r) calculated using Percus' test particle method. We 
find that there is reasonable qualitative agreement be- 
tween the equilibrium DFT density profiles displayed in 
Fig. Hm and the BD simulation results displayed in Fig. [21 
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FIG. 11: The free energy landscape F(w) plotted as a func- 
tion of w, the width of p a (r,w), calculated using both the 
DDFT and equilibrium DFT approaches, compared to the ex- 
act result, Eq. (|65[1 . For per 3 = 0.8 we find two disconnected 
branches of F(w). 

However, the profiles predict too much infilling in the re- 
gion close to the origin r < a, particularly at higher 
densities, which in turn results in an underestimate in 
the structure of the profiles at larger r > a. This error 
occurs for the reasons discussed in the previous subsec- 
tion IIV CI i.e. that the RY functional does not exert a 
strong enough interaction from the test particle onto the 
rest of the fluid. 

For per 3 = 0.8, shown in Fig. UPT d). we are able to 
calculate density profiles for all values of w, even though 
in the DDFT the profiles became 'trapped' at small w. 
The most striking aspect of these density profiles is that 
for intermediate values of w we see that p s (r) exhibits 
a plateau and a long tail. These features were not ob- 
served in the BD simulation results. However, similar 
features are present in G s (r, t) at intermediate times for 
colloidal spheres at densities close to the glass transition, 
where they are interpreted as the signature of dynamical 
heterogeneity in the system 3 -. 

D. Free energy landscape 

Since we are able to convert the dynamic density pro- 
files into their equilibrium equivalents via a set of external 
potentials, c.f. Eq. ((55)1. we are also able to calculate the 
equilibrium Helmholtz free energy for this correspond- 
ing equilibrium situation. Although this free energy is 
strictly an equilibrium quantity, since it underlies the 
time evolution of our dynamic approach we believe that 
it plays a relevant role. Therefore, by substituting the 
density profiles calculated using the DDFT into the free 
energy functional, Eqs. (|5Tj| and ((531) . we are able to map 
out a 'free energy landscape' as a function of t or w. Fig. 
[TT1 plots this free energy landscape, F(w), for densities 
per 3 = 0.2, 0.4, 0.6, 0.7 and 0.8. For per 3 = 0.2, 0.4, 
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FIG. 12: The free energy landscape F(w) plotted as a function 
of the width w/a, calculated by substituting the density pro- 
files from the exact equilibrium solution [Eqs. (|61l) and (|63[) ] 
into the RY functional. Here the deviation from the exact 
free energy and the emergence of the minimum are entirely 
due to the RY functional. 
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FIG. 13: The (scaled) strength of the parabolic potential, (3\, 
versus the resulting width, w/a, of p s (r). Results are calcu- 
lated using the equilibrium DFT and compared to the exact 
result. For the bulk fluid densities per 3 = 0.2, 0.4, 0.6 and 0.7, 
A decreases monotonically with w. However, for per 3 = 0.8 the 
curve is no longer monotonic, and exhibits two disconnected 
branches. The inset shows the effect of insufficient simula- 
tion run times on the equivalent situation studied using BD 
computer simulations, where t* is the simulation run time. 



and 0.6 we find that F(w) decreases monotonically with 
w. This decrease is initially steep and then the gradient 
begins to reduce as w increases. For pa 3 = 0.7, we find 
that after the initial steep descent the landscape devel- 
ops an almost constant plateau, but there is still a very 
small negative gradient. For per 3 = 0.8 the decrease is 
rapid and then the landscape terminates abruptly as the 
density profiles reach an arrested state. 

In the equilibrium case where Ud(r,w) — 0, one has 
an exact solution (see Appendix \X§ for the free energy 
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landscape as a function of the width, 

F(w) = F id HZ'n) ~ I In (jjjp) , (65) 

where is the ideal Helmholtz free energy, and Z' N is 
an irrelevant constant representing the partition function 
of the fluid when the test particle is located at the ori- 
gin. We plot Eq. (|6"5|) alongside the landscapes from the 
DDFT approach in Fig. [TTJ We find that this curve is 
located close to the DDFT landscape for per 3 = 0.2 but 
that the deviation grows with increasing density. 

We also calculate the free energy landscape via the 
equilibrium approach described in the previous subsec- 
tion and compare this to the results from both the 
dynamic approach and the exact equilibrium result in 
Fig. [TT] For the lower densities, per 3 = 0.2, 0.4 and 0.6, 
we find that there is good agreement between the DDFT 
and equilibrium DFT approaches. For per 3 = 0.7 we 
find that there is good agreement at low w, but around 
w/a = 0.7 there is a local maximum in F{w) in the equi- 
librium DFT results that is not present in the DDFT re- 
sults. For per 3 = 0.8 we find that the DDFT free energy 
landscape terminates abruptly at a fairly low value of w. 
However, since we can calculate the density profiles using 
our equilibrium approach for all widths we can therefore 
calculate F(w) for all w. We find a free energy landscape 
with two disconnected branches. Therefore, for a range 
of values of A we find two solutions to the Euler-Lagrange 
equations with different widths. Whether this is an in- 
dication of dynamic heterogeneity or an artefact of the 
functional is an interesting question. 

If we take the density profiles calculated using the ex- 
act route [Eqs. (|6"Tj) and (f6"3"|)] and substitute these into 
the RY functional, we find that the free energy curves, 
shown in Fig. [T^J do not follow the exact result (|6"5j) , but 
are very similar to those from the DDFT and DFT ap- 
proaches and even exhibit minima for per 3 = 0.7 and 0.8. 
Therefore, we must conclude that it is largely a property 
of the RY approximation for the free energy functional 
that generates these minima. 

In Fig. [13] we plot the exact relationship, Eq. (|62|) . be- 
tween the strength of the confining potential, A, against 
the width, w, and compare it to the results from the 
equilibrium DFT approach. At the lowest densities the 
equilibrium DFT approach closely follows the exact re- 
sult, but as the density is increased the width decreases 
for a given A. For per 3 = 0.8 we find that this curve is no 
longer monotonic and has two disconnected branches. 

Systems which are in a glassy or jammed state are 
by definition non-ergodic. A way of modelling non- 
ergodicity in Brownian dynamics is to measure the state 
of the system over too short a time frame. We demon- 
strate this effect by simulating a fluid where a single 
tagged particle is trapped in a parabolic potential well; 
c.f. Eq. (1BT)]) . If the radius of the well is sufficiently large 
and the simulation time is too small, then the particle is 
not able to fully explore the outer regions of the potential 
well. In the inset of Fig. [T3] we plot two results pertain- 



ing to this scenario where the fluid density is per 3 = 0.8, 
and t* is the simulation run time. The results for the 
longer simulation run time agree well with the exact re- 
sult, but results over the shorter time underestimate the 
width, particularly at the smaller values of A. A similar 
effect may exist in the DFT approach, where the non- 
ergodicity arises from not including the states where the 
particle is far from the origin. Recall that formally the 
density profile from DFT is an average over all possible 
states of the system. Using an approximate functional 
some of these states may be neglected^. 



VI. CONCLUDING REMARKS 

On the basis of the theoretical and simulation results 
that we presented in this paper for the dynamics of the 
bulk hard sphere fluid, we conclude that the dynamical 
test particle limit, combined with DDFT, provides a re- 
liable method for calculating the van Hove (and other 
related) dynamical pair correlation functions at low and 
at intermediate densities per 3 < 0.6. In the previous 
publication^ we have shown that the theory may be 
applied in a fairly straightforward manner also to inho- 
mogeneous situations, hence we conclude that the dy- 
namic test particle theory may indeed be used to calcu- 
late the van Hove function for fluids at interfaces and un- 
der confinement. Furthermore we have shown that, quite 
surprisingly, at low and intermediate densities the very 
simple Vineyard approximation (reviewed in Sec. Ill Cj) is 
actually quantitatively fairly good. This approximation 
only requires as input the radial distribution function 
g(r) and the diffusivity Di and therefore provides a very 
useful quick approach for obtaining an approximation for 
the van Hove function for colloidal fluids. 

The possible conclusions about the performance and 
even about the qualitative status of the predictions of 
the theory at higher densities are far more intricate 
though. In this regime, the theory in its current form 
is clearly quantitatively unreliable - compare for exam- 
ple, the results from BD simulations in Fig. [2jd) to 
those from the dynamic test particle theory shown in 
Fig. [3Jd) for per 3 = 0.8, where the theory predicts that 
the system jams, whereas in simulations the system is 
an ergodic fluid. Moreover, at even higher densities, the 
monodisperse hard sphere system crystallises in simula- 
tions rather than undergoing a glass transition; polydis- 
persity would be required to suppress freezing^. We be- 
lieve that the behaviour of the theory at higher densities 
can primarily be ascribed to our choice of approxima- 
tion for the free energy functional - see Sec. lIVBl for the 
reasons for using the RY functional ((55)) in the present 
study. We believe that if we had used a more accu- 
rate functional, such as Rosenfeld's fundamental measure 
theory^ - — , quantitatively more accurate results could 
be obtained at these higher densities. Nevertheless, the 
theory yields a clear prediction that there is a dynamic 
(glass) transition, whereby the tagged (self) particle be- 
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comes trapped in the cage formed by the surrounding 
particles. We can offer three possible interpretations of 
this result. 

Firstly, one may conclude that this glass transition 
stems simply from the use of the approximate RY func- 
tional and since the theory predicts the glass transition 
to be at a density value that is well below where the 
glass transition is believed to occur, our results at higher 
densities should be disregarded. 

An alternative conclusion that one may draw is that 
the theory is correctly describing (some of) the physics 
of the glass transition, but that the predictions are only 
qualitative in nature and occur in reality at higher densi- 
ties and in polydisperse systems. Support for this point 
of view comes in particular from results such as those 
in Fig. [51 where we display results for 0(fcer = 2ir,t). 
For pa 3 = 0.7, which is near to where the theory pre- 
dicts the glass transition to be, we find a plateau in 
4>(k, t) and the clear presence of two-stage relaxation, 
which indeed has been observed in hard sphere colloidal 
suspensions 9 . Hence our results are (qualitatively) simi- 
lar to those from MCT-' 9 . Further support for the above 
interpretation comes from the results in Figs. [TT1and[T21 
for the behaviour of the free energy landscape that un- 
derlies the DDFT. The appearance of a minimum in the 
free energy as a function of displacement, correspond- 
ing to a particle becoming trapped in the cage formed 
by its neighbours, is one central prediction of the the- 
ory by Schweizer and Saltzmann— who combined ele- 
ments of MCT, DFT, and activated rate theory, in order 
to describe localization and transport in glassy fluids. 
Furthermore, given some of the work in the literature 
based on DFT to study the glass transition, our predic- 
tion that particles become localized should not come as 
a surprise: Wolynes and coworkers^— developed a suc- 
cessful model of hard sphere vitrification, which is similar 
to the DFT treatment of crystallization^^. Using a ran- 
dom close-packed, non-periodic lattice they found a fluid- 
glass transition where the fluid "crystallizes" onto this 
lattice. The success of this method, along with its ability 
to model the freezing transition (onto a regular lattice), 
has provoked a number of further developments^ 7 - - — . 
Other approaches^ have investigated dense Brownian 
systems through modelling via certain stochastic differen- 
tial equations and found that the system exhibits glassy 
behaviour. Thus, overall our results seem to be qualita- 
tively consistent with other DFT based theories and with 
MCT for the glass transition. Nevertheless, the density 
where the glass transition occurs, as predicted by the the- 
ory in its present form, is far too low. Furthermore, it 
could be the case that the similarity between our results 
and those from the MCT arc somehow a mathematical 
(rather than physical) coincidence, since an essential fea- 
ture of MCT is the presence of memory in the dynamical 
equations. This important feature is absent from the 
present DDFT. 

The third possible conclusion that one may draw con- 
cerns the question whether the minimum in the free en- 



ergy and the localization of the tagged self particle are 
merely a signature of freezing in the theory. The RY func- 
tional is well-known to predict the freezing transition to 
occur at a density below that where it occurs in reality. 
It could simply be the case that this functional overly 
favours freezing, so that when it is applied in the way 
we use it here, where we constrain all density profiles, 
p s (r) and Prf(r), to be spherically symmetric, a signature 
of freezing shows up as the tagged self particle becoming 
localized. 

Some merit can be found in all of the arguments out- 
lined above and we find ourselves unable to judge which 
one(s) are correct. Indeed further work is required to 
provide a clear assessment of these issues. In particu- 
lar, the dynamical test particle theory should be imple- 
mented with a more sophisticated approximation for the 
free energy functional than we have used here. 

As we have shown, our approach is based on integrat- 
ing the Smoluchowski equation ([^)) over all except one 
of the position coordinates, in order to derive an equa- 
tion for the one-body density distribution. An alterna- 
tive approach is to integrate over all but two of the po- 
sition coordinates, in order to obtain (|23p an equation 
for the two-body distribution function p^ (ri,r2,i) = 
N(N — 1)J fdr 3 ...fdr N P(r N ,t). The resulting dy- 
namical equation depends on the three body distribution 
function p( 3 )(ri, r 2 , r 3 , t). On making a suitable closure 
approximation, this provides a different starting point for 
studying the pair correlations in a colloidal fluid - see e.g. 
Refs. [52II53I and references therein, which also consider 
the effect of the hydrodynamic interactions between the 
colloids. Developing the theory for the dynamical pair 
correlation functions in this way is very natural. How- 
ever, we believe that the strength of our method, where 
we use the dynamical test particle approach allowing us 
to work at the one-body level, is that we are able to use 
DFT to close our equations and therefore we are able to 
describe the fluid spatial correlations very accurately. 

Finally, we mention other possible directions for de- 
veloping the theory in the future. One important aspect 
in the dynamics of colloidal dispersions, that we have 
entirely neglected here, are the hydrodynamic interac- 
tions between the particles. Rex and Lowe n 73 ' 74 have 
shown how to include the hydrodynamic interactions in 
a DDFT treatment and so it would be worthwhile to use 
their DDFT formulation together with the present dy- 
namical test particle limit, in order to calculate the van 
Hove function under the influence of hydrodynamic in- 
teractions. 

A further aspect of our work, that offers possible exten- 
sions of the theory, concerns the question how to model 
the diffusivity of the tagged particle in a better way: In 
the dynamical test particle calculation one could replace 
the (constant) diffusion coefficient in Eq. (J3FJ with a dif- 
fusion coefficient that depends on time; i.e. to replace 
D — > D (t). In doing this one could ensure that D(t) 
takes the correct values at both short and long times. 
However, doing this still does not treat memory effects 
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in the dynamics. As MCT demonstrates, memory effects 
are key for a system to exhibit the ideal glass transition 
Thus, we believe that including memory 
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scenari* 

into our theory would be a crucial step in future work. 
This could possibly be done along the lines of the in- 
teresting work of Mcdina-Noyola and coworkers^— . To 

include memory in our theory one could replace Eq. (29) 
withi^S: 



dp l {r,t) 
dt 



V • J dt' J drT(r-r',t-t') 

SF[{piY 



6 Pi (r>,t>) 



(66) 



where the mobility coefficient T has been replaced by 
one that is non-local in time and space. However, this 
would result in a considerable increase in computational 
complexity as within DFT the correlations in space are 
already treated in a complex manner and these would 
need to be coupled to the correlations in time. Whether 
such non-locality helps to cure some of the deficiencies of 
our approach is an open question. 



our external potential (|A3I) into (|A5I) we obtain 

Zn = y"dr ie xp(-/3Ar2) J dr 2 .. N exp[-j3V(r N )) 
= J dr lG xp(-/3Ar 2 ) J dv' 2 .. N exp(-pV(r' 2 .. N )), 

where in the second step we have made the substitution 
r- = rj — ri, for i = 2..N, so that we can do the integra- 
tions over the positions r' 2 ...r' N . This gives, 

Z N =/dr 1 exp(-/3Ar?)Zj v = (tt / (3\f 2 Z' N , 

where Z' N is the configuration integral for N particles 
where one particle is located at the origin. The Helmholtz 
free energy is then given by, F = — ji~ x ]n(Qjv(V, T)) = 

— ln(Qjy Zn v^ T%> ) where Vjv is the volume occupied 
by the particles, which yields 



PF = F id -ln(Z' N /V N )-^\n(jj- 



(A6) 



Appendix A: Exact Results 

We consider a fluid of N particles with positions rj, 
momenta Pi and mass m in the presence of an arbitrary 
external field that acts only on particle i = 1, ui(r) = 
Ar 2 . Assuming that we are in the classical limit, the 
Hamiltonian is given by Hn = K + V + U where the 
contributions are due to the (classical) kinetic energy, 
the total inter-particle potential (not necessarily pairwisc 
additive), and the external potential, respectively; 



N 

k = 



2m 



V = u(r 1( ..,rjv) 



(Al) 

(A2) 
(A3) 



U = \v{ 

The canonical partition function, (5jv(V,T), is given by 
h' 3N 



Qn{V,T) 



(N-iy. 



dr^dp- eM-PHNir",^)} 
(A4) 

where h is Planck's constant, and the (N — 1)! factor 
results from the fact that besides particle i = 1, the re- 
maining particles are indistinguishable. The integrations 
over momenta in Eq. (|A4I) can be carried out explicitly, 
leaving a configuration integral over positional degrees of 
freedom: 



/ 



dr 1 ...dr Ar exp(-/3(F + C/)). (A5) 



Note that for Brownian particles Zn is also the quantity 
that characterises the structure of the fluid. Substituting 



Therefore, the Helmholtz free energy only depends on the 
confining potential in a simple way. 

One can also obtain the one body density profiles. In 
general, for a system of N particles, the one body density 
profile, pjy (r) can be obtained fromi, 



ZNi N-l)l J CXP [-^(r^) + *(r w ))] , 

(A7) 



where the Nl/(N — 1)! factor accounts for the indistin- 
guishability of the particles. 

For the single particle subject to the external potential 
we get 



Pi 1] (ri) 



exp(— /3Arf ) 



dr 2 .. Jv exp(~/3F(r A ')), 
dr^expH3V(r^jv)), 



exp(— f3Xr\ ) 
(n/f3X)^Z' N 
exp(— /3Arf ) 

(7T//3A) 3 / 2 ' 



which is a normalised Gaussian. It can be shown that 
since p s (r) is a Gaussian, then w and A are simply related 
by 



w = (2A/3)- 1 / 2 , 
and we can rewrite Eq. (| A6|) as 



(A8) 



(A9) 



We now seek the density profile of the remaining parti- 
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cles: 



4 1} (r 2 ) 



x J dr 3 .. N cM-PV(r N )), 

J dri exp(— /3Ar 2 ) 
(tt//?A) 3 / 2 
(A-l) 



(A10) 



dra.^expC-^^)) 



To progress we make use of the formal relationship be- 
tween g[r) and the two body density profile, pw (ri,^), 
which for a homogeneous fluid can be shown to bei, 



SfV(ri-r 2 ) = -pPk {ri-r 2 ), 
1 Nl 



p 2 Z N (N-2)\ 

x y dr 3 .. Jv exp(-/31/(r A ')), 
V N(N- 1) 

x J dr 3 .. Jv exp(-/3V(r JV )), 



where we have made the substitutions, p — N/V and 
Zn — Z' N V. Cancelling terms and rearranging we get, 



p<4 2) (n - r 2 ) = { -?LJ± I dr 3 .. N exp(-pV(r N ))- 



Z' 



N 



(All) 



Substituting (|Allj) into (|A10|) gives, 

(1) Jdriexp(-^Ar 2 ) (2) 

« (r2) = (7r// 3A)3/2 PgV(n-r 2 ) 

= (^5372 / dneipH^a^n-is), 

which is the normalised Gaussian convolved with pg(r). 
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